function B=MagneticField3D(x_d,y_d,z_d,l,ai,ao)

r_d=sqrt(x_d.^2+y_d.^2);

if r_d>0
    cos_theta_d=x_d/r_d;
    sin_theta_d=y_d/r_d;
    
    B_R=triplequad(@(r,z,theta)MagneticField3D_R_Func...
        (r,z,theta,r_d,z_d),ai,ao,-l/2,l/2,0,2*pi);
else
    B_R=0;
    cos_theta_d=0;
    sin_theta_d=0;
end

B_Z=triplequad(@(r,z,theta)MagneticField3D_Z_Func...
    (r,z,theta,r_d,z_d),ai,ao,-l/2,l/2,0,2*pi);

B=-[B_R*cos_theta_d B_R*sin_theta_d B_Z];
